Периодограмма

Предположим, что временной ряд образован линейной комбинацией синусоид и косинусоид взятых с различными частотами. Предложил эту модель был Шустер (1898 год) В это время очень был популярен анализ Фурье и Шустер был первым, кто применил методы рядов Фурье к исследованию стохастических временных рядов. Для начала предположим, что длина ряда N нечетна N= 2q+1. Методом наименьших квадратов будем подгонять к исходному ряду модель ряда Фурье.

\(x_{t}=\alpha_0+\sum_{i=1}^{q}(\alpha_{i}c_{it}+\beta_{i}s_{it})+\epsilon_{t}\)

Где

\(c_{it}=cos2π\lambda_{i}t , s_{it}=sin2π\lambda_{i}t, \lambda_{i}=(i/N)\)

Нетрудно убедиться, что оценками неизвестных параметров \(\lambda_{i}\) и \(\beta_{i}\) будут

\(\alpha_0=\sum_{t=1}^{N}x_{i}=x\)

\(\alpha_{i}=(2/N)\sum_{t=1}^{N}x_{t}c_{it}, i=1,...,q\)

\(\beta_{i}=(2/N)\sum_{t=1}^{N}x_{t}s_{it}, i=1,...,q\)

Рассмотрим функцию вида

\(I(\lambda_{i})=(N/2)(a_{i}^2+b_{i}^2), i=1,...,q\)

Эта функция получила название периодограммы. В случае, если \(N\) четно и \(N= 2q\), то приведенные выше формулы для \(i=1,...,q-1\) не изменятся, а при \(i=q\) примут вид

\(\alpha_{q}=(1/N)\sum_{t=1}^{N}(-1)^{t}x_{t}\)

\(\beta_{q}=0\)

Для примера рассмотрим детерминированный ряд

\(y_t= 2cos(2\pi t\frac{4}{96})+3cos(2\pi t\frac{14}{96}+0.3)\)

t<- 1:96
cos1<- cos(1*pi*t*4/96)
cos2<- cos(1*pi*t*14/96+0.3)
y <-2*cos1+3*cos2
plot(t,y, type = "o",ylab= expression(y[t]),col = "blue",lwd = 2)

Периодограмма для ряда \(y_t\)

library(TSA)

## Loading required package: leaps

## Loading required package: locfit

## locfit 1.5-9.1    2013-03-22

## Loading required package: mgcv

## Loading required package: nlme

## This is mgcv 1.8-9. For overview type 'help("mgcv-package")'.

## Loading required package: tseries

## 
## Attaching package: 'TSA'

## The following objects are masked from 'package:stats':
## 
##     acf, arima

## The following object is masked from 'package:utils':
## 
##     tar

periodogram(y,col = "blue",lwd = 2)
abline(h=0)

Теперь добавим шум. Для двух случайно выбранных частот из набора частот \([1/96,2/96,...,47/96\)], амплитуды \(a,b\) также выберем случайно и добавим белый шум Итак модель будет

\(y_t= A_1cos(2\pi f_1 t)+B_1sin(2\pi f_1t )+A_2cos(2\pi f_2 t)+B_2sin(2\pi f_2t )+\epsilon_t\)

set.seed(134); t=1:96; integer=sample(48,2)
freq1=integer[1]/96; freq2=integer[2]/96
A1=rnorm(1,0,2); B1=rnorm(1,0,2)
A2=rnorm(1,0,3); B2=rnorm(1,0,3); w=2*pi*t
y=A1*cos(w*freq1)+B1*sin(w*freq1)+A2*cos(w*freq2)+
B2*sin(w*freq2)+rnorm(96,0,1)
plot(t,y,type='o',ylab=expression(y[t]),col="blue",lwd=2)

Периодичность и частоты этих периодичностей весьма не очевидны. Построим периодограмму

periodogram(y, col = "blue",lwd = 3)
abline(h=0)

Где выбранные частоты становится очевидно.

Cпектральная плотность

При определении периодограммы мы предполагали, что частоты \(\lambda_{i}=(i/N)\). При переходе к спектру мы ослабим это предположение, считая, что частота \(\lambda\) непрерывно в диапазоне от 0 до 0.5. При этом определение периодограммы примет вид

\(I(λ)=(N/2)(a_{\lambda}^2+b_{\lambda}^2), \lambda∈[0,0.5]\)

Эту функцию будем называть выборочным спектром. Справедливо следующее очень важное утверждение.

Теорема 1. Выборочный спектр - косинус-преобразование Фурье автоковариационной функции. Доказательство

Обозначим через \(d _{ \lambda}\) - комплексное число равное разности оценки амплитуды косинусоидальной и синуcоидальной составляющей умноженной на \(i=\sqrt{-1}\)

\(d _{ \lambda}=a_{ \lambda}-ib_{ \lambda}\)

Тогда выражение для выборочного спектра можно переписать следующим образом.

\(I(\lambda)=(N/2)(a_{\lambda}^2+b_{\lambda}^2)=(N/2)(a_{\lambda}-ib_{\lambda})(a_{\lambda}+ib_{\lambda})=(N/2)d_{\lambda}d_{\lambda}'\)

Из формул для оценок наименьших квадратов имеем

\(d_{\lambda}=(2/N)[\sum_{t=1}^{N} x_{t}c_{\lambda}t -i\sum_{t=1}^{N}x_{t} s_{\lambda}t]=(2/N)[\sum_{t=1}^{N} x_{t}cos(2\pi\lambda t) -i\sum_{t=1}^{N}x_{t} sin(2\pi\lambda t]= (2/N)\sum_{t=1}^{N} x_{t}e^{-2\pi\ i\lambda t}=\)

\((2/N)\sum_{t=1}^{N} (x_{t}-\overline x)e^{-2\pi\ i\lambda t}\)

Подставляя полученное выражение в формулу для спектра, получим

\(I(\lambda)=(N/2)d _{\lambda}d _{\lambda}'= (N/2)( (2/N) \sum_{t=1}^{N} (x_{t}-\overline x)e^{-2\pi i \lambda t}) ((2/N)\sum_{t'=1}^{N}(x_{t}-\overline x)e^{-2\pi i \lambda t′})=\)

\((2/N)\sum_{t=1}^{N}\sum_{t′=1}^{N}(x_{t}-x)(x_{t′}-x)e^{-2\pi i\lambda (t-t′)}\)=

Меняя порядок суммирования, и ,обозначив \(k=t-t′\), последнее выражение можно переписать

\(2/N\sum_{k=-N+1}^{N+1}e^{-2\pi i\lambda k} \sum_{t=1}^{N-k}(x_t-\overline x)(x_{t+k}-\overline x)\)

Так как

\(c_k =\frac{\sum_{t=1}^{N-k}(x_t-\overline x)(x_{t+k}-\overline x)}{N}\)

Получим

\(I_{N}(\lambda)=2\sum_{k=-N+1}^{N-1}с_{k}e^{-2\pi i\lambda k}\) \(=2[c_0+2\sum_{k=1}^{N-1}c_{k}cos(2\pi\lambda k)]\)

Спектральное представление автоковариационной функции случайного процесса

Пусть \(x_t\) стационарный случайный процесс, \(c_k, k = 0,\pm 1,\pm 2,...,\pm \infty\) его автоковариационная функция. Справедливо следующее представление для АКФ стационарного случайного процесса

\(с_k=\int_{-0.5}^{0.5}e^{2\pi\lambda k}dF(\lambda)\)

Это выражение называется спектральным представлением АКФ. Функция \(F(\lambda)\) называется спектральной функцией случайного процесса. В случае абсолютной непрерывности \(F(\lambda)\) последнее выражение записывают как

\(с_k=\int_{-0.5}^{0.5}e^{2\pi i\lambda k}f(\lambda)d\lambda\)

где \(f(\lambda)\)- называют спектральной плотностью.

И наоборот

\(f\lambda) = \sum_{k=-\infty}^{\infty}c_k e^{-2\pi i\lambda k}= c_0+2\sum_{k=1}^{\infty}c_kcos(2\pi\lambda k)\)

Свойства спектра

\(F(\lambda)\) не убывает
непрервывна справа

3.\(F(\lambda)\) >0

\(\lim_{\lambda -> 1/2} F(\lambda)= D[(x_t)] = c_0\)

Спектральная плотность ARMA процесса

Спектральная плотность линейного фильтра

Пусть \(x_t\) некоторый cтационарный временной ряд и \(...,c_{-1},c_0,c_1,c_2,...\) бесконечная абсолютно суммируемая последовательность чисел. Построим новый временнной ряд по правилу

\(y_t= \sum_{k=-\infty}^{\infty}c_kx_{t-k}\)

это модель линейного фильтра, это выражение называют также дискретной сверткой последовательностей \({c_k}\) и \({x_t}\). Пусть \(f_x(\lambda),f_y(\lambda)\) спектральные плотности процессов \(x_t, y_t\) соответственно. Пусть

\(C(e^{-2\pi i \lambda})= \sum_{j=-\infty}^{\infty}c_j e^{-2\pi i\lambda j}\)

Тогда

\(cov(y_t,y_{t-k})= cov(\sum_{j=-\infty}^{\infty}c_jx_{t-j},\sum_{s=-\infty}^{\infty}c_sx_{t-k-s})=\)

\(\sum_{j=-\infty}^{\infty}\sum_{s=-\infty}^{\infty}c_jc_scov(x_{t-j},x_{t-k-s)}\)

но

\(cov(x_t,x_{t-k})=\int_{-0.5}^{0.5}e^{2\pi i\lambda k}f_x(\lambda)d\lambda\)

поэтому

\(cov(y_t,y_{t-k})= \sum_{j=-\infty}^{\infty}\sum_{s=-\infty}^{\infty}c_jc_s\int_{-0.5}^{0.5}e^{2\pi i\lambda (s+k-j)}f_x(\lambda)d\lambda=\) \(\int_{-0.5}^{0.5}|\sum_{s=-\infty}^{\infty}c_se^{-2\pi i\lambda s}|^2 e^{2\pi i \lambda k}f_x(\lambda)d\lambda\)

Итак

\(cov(y_t,y_{t-k})=\int_{-0.5}^{0.5}|C(e^{-2\pi i\lambda s})|^2 e^{2\pi i \lambda k}f_x(\lambda)d\lambda\)

но по определению

\(cov(y_t,y_{t-k})= \int_{-0.5}^{0.5}e^{2\pi i\lambda k}f_y(\lambda)d\lambda\)

В результате получим

\(f_y(\lambda)= |C(e^{-2\pi i\lambda})|^2f_x(\lambda)\)

Функцию \(|C(e^{-2\pi i\lambda})|^2\) называют переходной функцией фильтра

Спектральная плотность белого шума

Непосредственно из выражения

\(f\lambda) = c_0+2\sum_{k=1}^{\infty}c_kcos(2\pi\lambda k)\)

следует, что для белого шума \(f(\lambda)= c_0=\sigma_{\epsilon}^2\) это объясняет слово “белый” в названии шума.

Спектральная плотность MA(1) процесса

MA(1) процесс это модель линейного фильтра с \(с_0=1\) и \(с_1=-\theta\) Поэтому

\(|C(e^{-2\pi i\lambda})|^2= (1-\theta e^{2\pi i\lambda})(1-\theta e^{-2\pi i\lambda})= 1+\theta^2-\theta (e^{2\pi i\lambda}+e^{-2\pi i\lambda})= 1+\theta^2-2\theta cos(2\pi\lambda)\)

Поэтому для MA(1) процесса спектральная плотность выглядит следующим образом

\(f(\lambda)= [1+\theta^2-2\theta cos(2\pi\lambda)]\sigma_{\epsilon}^2\)

theta = 0.9
ARMAspec(model = list(ma = -theta))

Для произвольных \(\theta\) вид спектральной плотности можно посмотреть в следующем фрейме.

Спектральная плотность AR(1) процесса

Модель авторегрессии первого порядка это тоже модель линейного фильтра только на вход идет AR(1) процесс а на выход белый шум, поэтому выражение для спектральной плотности

\(f(\lambda)[1+\theta^2-2\theta cos(2\pi\lambda)]= \sigma_{\epsilon}^2\)

или

\(f(\lambda)= \frac{\sigma_{\epsilon}^2}{[1+\phi^2-2\phi cos(2\pi\lambda)]}\)

phi = 0.9
ARMAspec(model = list(ar = phi))

Для произвольных \(\phi\) вид спектральной плотности можно посмотреть в следующем фрейме.

Спектральная плотность AR(2) процесса

Аналогично можно получить вид спектральной плотности для процесса AR(2) \(f(\lambda)= \frac{\sigma_{\epsilon}^2}{[1+\phi_1^2+\phi_2^2- 2\phi_1(1-\phi_2)cos(2\pi\lambda)- 2\phi_2cos(4\pi\lambda)]}\)

phi_1 = 0.5
phi_2 = -0.4

ARMAspec(model = list(ar = c(phi_1,phi_2)))

Спектральная плотность ARMA(1,1) процесса

Для смешаных процессов вид спектральной плотности

\(f(\lambda)= \frac{\sigma_{\epsilon}^2[1+\theta^2-2\theta cos(2\pi \lambda)]}{[1+\phi^2-2\phi cos(2\pi\lambda)]}\)

Подробнее с поведением спектральной плотности смешаного процесса можно ознакомиться в следующем фрейме.

Спектральная плотность ARMA(p,q) процесса

Для смешаных процессов произвольно порядка спектральную плотность удобно записывать через характеристические полиномы

\(S(\lambda)= |\frac{\theta_q(e^{-2\pi i \lambda})}{\phi_p(e^{-2\pi i \lambda})}|^2\sigma_{\epsilon}^2\)

phi_1 = 0.5
phi_2 = -0.4
theta_1 = -0.3
theta_2 = 0.6
ARMAspec(model = list(ar = c(phi_1,phi_2),ma = c(theta_1,theta_2)))

Ceзонные процессы

Здесь также будем использовать соответствующие характеристические полиномы \(\phi_p(B),\theta_q(B),\Phi_P(B),\Theta_Q(B)\). Спектральная плотность сезонного процесса задается выражением

\(f(\lambda)=\sigma_{\epsilon}^2|\frac{\Theta_Q(e^{-2\pi i \lambda})\theta_q(e^{-2\pi i \lambda})}{\Phi_P(e^{-2\pi i \lambda})\phi_p(e^{-2\pi i \lambda})}|^2\)

Для примера посмотрим поведение спектральной плотности сезонного процесса авторегрессии

\((\Phi_1 B^S+\Phi_2 B^{2S})(\phi_1 B+\phi_2 B^2)x_t=\epsilon_t\)

Оценка спектральной плотности по выборке

В качестве введения в проблемы, возникающие при оценке спектральной плотности, рассмотрим следующий численный пример. Пусть смоделирован процесс AR(1) c параметром \(\phi=-0.7\).Длина выборки равна 200. Нам известны его свойства и вид спектральной плотности.

n <- 200
phi <- -0.7
y <-arima.sim(model = list(ar = phi),n=n)
sp <- spec(y,col = "blue",lwd = 2,xlab= "frequency",ylab = "spectral density")
lines(sp$freq,ARMAspec(model = list(ar = phi),freq = sp$freq,plot= "F")$spec,lty = 2)
abline (h = 0)

Для изучения свойств оценки спектральной плотности предположим, что \(y_t\) есть нормальный белый шум с дисперсией \(c_0\). Мы уже знаем выражения для оценки амплитуд для частот \(\lambda = j/n\)

\(A_{\lambda}= 2/n\sum_{t=1}^n y_t cos(2\pi t \lambda)\)

\(B_{\lambda}= 2/n\sum_{t=1}^n y_t sin(2\pi t \lambda)\)

\(A_{\lambda}\) и \(B_{\lambda}\)- линейные функции от \(y_t\) следовательно обе имеют нормальное распределение. Можно оценить их математическое ожидание и дисперсию. Среднее будет 0, дисперсия \(2c_0/n\). Более того так синус и косинус ортогональны то \(A_{\lambda}\) и \(B_{\lambda}\) будут некоррелированы, а значит в силу нормальности независимы для каждого \(\lambda\). Сложнее показать, но поверим, что взаимно независимы будут \(A_{\lambda_1},A_{\lambda_2},...\) и \(B_{\lambda_1},B_{\lambda_2},...\). Квадрат нормальной сл. величины имеет хи-квадрат распределение с одной степенью свободы. Сумма независимых хи-квадрат имеет также хи-квадрат распределение с числом степеней свободы равной сумме степеней свободы каждой случайной величины. Таким образом

\(n/2[A_{\lambda}^2+B_{\lambda}^2]= \frac{2\overline I(\lambda)}{f(\lambda)}\)

имеет хи-квадрат распределение с двумя степенями свободы. Также имеем , что \(\overline I_{\lambda_1}\) и \(\overline I_{\lambda_2}\) независимы при \(\lambda_1 \ne \lambda_2\). И

\(E[\overline I_{\lambda}]= f({\lambda})\)

\(D[\overline I_{\lambda}]= f({\lambda})^2\)

К сожалению последние выражения не зависят от длины выборки \(n\). Это означает, что выборочная спектральная плотность не являетя состоятельной оценкой теоретической спектральной плотности. Увеличение длины выборки не приведет к более гладкой оценке. Одним из возможных решений, но неудобных на практике, это повторить моделирование много раз и усреднить оценку

n <- 200
phi <- -0.7
y <-arima.sim(model = list(ar = phi),n=n)
sp <- spec(y,col = "blue",lwd = 2,xlab= "frequency",plot = FALSE,ylab = "spectral density")
totalspec <- sp$spec
m <- 40
for (i in 2:m)
{
  y <-arima.sim(model = list(ar = phi),n=n)
  sp <- spec(y,col = "blue",lwd = 2,xlab= "frequency",plot = FALSE,ylab = "spectral density")
  totalspec <- totalspec + sp$spec
}
totalspec <- totalspec/m
plot(sp$freq,totalspec,type = "l", col = "blue", lwd = 3)
lines(sp$freq,ARMAspec(model = list(ar = phi),freq = sp$freq,plot= "F")$spec,lty = 2,lwd =2,col = "red")
abline (h = 0)

Картина становится лучше,но на практике такой прием мало применим. Чаше всего имеется одна выборка и необходимо получить достаточно гладкую оценку плотности.

Сглаженная оценка спектральной плотности

Основная идея ,так как спектральная плотности не сильно меняется при малом изменении частоты, то мы будем усреднять оценку на близких частотах. Итак мы пришли к понятию сглаженной спектральной плотности. Зададим некоторое произвольное число \(m>0\), которое будем называть шириной окна и с весами \(\frac{1}{2m+1}\) усредним выборочную спектральную плотность

\(\overline{f}(\lambda)= \frac{1}{2m+1} \sum_{j=-m}^mf(\lambda+j/n)\)

Другими словами мы сгладили спектральную плотность с помощью спектрального окна \(W_m(\lambda)\) ширины \(m\) свойства которого.

\(1.W_m(k)>0\)

\(2.W_m(k)=W_m(-k)\)

\(3.\sum_{k=-m}^m W_m(k)=1\)

и получили сглаженную спектральную плотность

\(\overline{f}(\lambda)= \sum_{j=-m}^m W_m(k)f(\lambda+j/n)\)

Веса \(\frac{1}{2m+1}\) называют простым спектральным окном, еще его называют спектральным окном Даниелла по имени того, кто первым использовал это в 1940 году. Для примера рассмотрим сглаженное спектральное окно Даниелла ширины 5 для смоделированного AR(1) процесса

library(TSA)
set.seed(271435); n=200; phi=-0.6
y=arima.sim(model=list(ar=phi),n=n)
k=kernel('daniell',m=5)
sp=spec(y,kernel=k,log='no',sub='',xlab='Frequency',
ylab='Smoothed Sample Spectral Density',col = "blue",lwd = 2)
lines(sp$freq,ARMAspec(model=list(ar=phi),freq=sp$freq,
plot=F)$spec,lty='dotted')

\(D[\overline{f}(\lambda)]=\sum_{k=-m}^m W^2_m(k)D[f(\lambda+k/n)]=\sum_{k=-m}^m W^2_m(k)f^2(\lambda)\)

Таким образом

\(D[\overline{f}(\lambda )]=f^2(\lambda) \sum_{k=-m}^m W^2_m(k)\)

В частности для окна Даниелла

\(\sum_{k=-m}^m W^2_m(k)= \frac{1}{2m+1} -> 0\) при \(n->\infty\)и \(m->\infty\) и сглаженная оценка спектральной плотности cтановится состоятельной.

В последующие годы было предложено масса других спектральных окон. Один из подходов здесь неоднократное повторение процедуры сглаживания с одним и тем же окном в частности Даниелла. Это как бы свертка окон.

k<- kernel('daniell',5)
k2<- kernel('daniell',c(5,5))
plot(k,type= "h",lwd = 4,col = "blue")

plot(k2,type= "h",lwd = 4,col = "blue")

k3<- kernel('daniell',c(5,5,5))
plot(k3,type= "h",lwd = 4,col = "blue")

Сгладим ранее смоделированный ряд с этим окном и сравним с прошлыми результатами сглаживания и выборочной периодограммой

sp3=spec(y,kernel=k3,log='no',sub='',xlab='Frequency',
ylab='Smoothed Sample Spectral Density',col = "blue",lwd = 2)
lines(sp$freq,ARMAspec(model=list(ar=phi),freq=sp$freq,
plot=F)$spec,lty='dotted')
lines(sp$freq,sp$spec,lwd = 2,col = "red")
legend("topleft",c("Daniell window","Daniell^3 window","True spectrum","Periodogram"),bty="n",lwd = 2,col = c("red","blue","black","magenta"))
per <- periodogram(y,plot = FALSE)
lines(per$freq,per$spec,col = "magenta")

Выбор ширины окна важная задача. Чатфилд (Chatfield) 2004 предложил выбирать ширину окна \(m = \sqrt{n}\) Другие предлагают вариант \(m=2\sqrt{n}\)

Корреляционное окно

Здесь мы вспомним, что спектр это косинус преобразование Фурье автокорреляционной функции

\(\overline{f}(\lambda)= \sum_{k=-m}^m W(k)f(\lambda+k/n)=\)

\(\sum_{k=-m}^m W(k)[\sum_{j=-n+1}^{n-1}c_j e^{-2\pi i (\lambda+k/n)j}]=\)

\(\sum_{j=-n+1}^{n-1}c_j[\sum_{k=-m}^m W(k)e^{-2\pi i k/nj}]e^{-2\pi i \lambda j}\)

Выражение в квадратных скобках обозначим через \(w(j/n)\)

\(w(j/n)= \sum_{k=-m}^m W(k)e^{-2\pi i k/nj}\)

Итак

\(\overline{f}(\lambda) = \sum_{j=-n+1}^{n-1}w(j/n)c_j cos(2\pi\lambda j)\)

функция \(w(x)\) имеет свойство

\(w(x)= w(-x)\)

\(w(0)= 1\)

\(w(x) <= 1\)

поэтому \(w(x)\) называют корреляционным окном

Примеры корреляционных окон

1.Прямоугольное окно \(w(j/n)= 1\) при \(|j| <= m\) Здесь \(m\) ширина корреляционного окна

Треугольное окно

\(w(j/n)= 1- |j/m|\) при при \(|j| <= m\)

Энцефалограмма головного мозга (ЭЭГ)

Наблюдения снимаются аппаратом ЭКГ 128 раз в секунду, т.е. интервал между наблюдением равен 0.0078 сек.

eeg <- read.csv("c:/ShurD/aLection/EEG.csv")
matplot(eeg,type ="l",col = "blue",main= "EEG")

k3<- kernel('daniell',c(5,5,5))
eegspec=spectrum(eeg,kernel=k3,log='no',sub='',xlab='Frequency',main = "EEG.Periodogram and Smoothed Spectral Density",
ylab='Spectral Density',col = "blue",lwd = 2)
pereeg <- periodogram(eeg,plot = FALSE)
lines(pereeg$freq,pereeg$spec,col = "magenta")  
legend("topright",c("Daniell^3 window","Periodogram"),bty="n",lwd = 2,col = c("blue","magenta"))

Первый пик слева находится на частоте 0.02 второй на 0.09, что соответсвует периодам 1/0.02 = 50, второй 1/0.09= 11.1. В секундах периоды равны соответственно 50 * 0.0078 сек.= 0.39 сек. и 11.1 *0.0078 = 0.086 сек.. Выразим их в герцах 1/0.39 = 2.7 Гц. и 1/0.86 = 11.5 Гц. \(\alpha\) ритм мозга имеет частоты 9-13 Гц. ;\(\delta\) ритм 1-3 Гц.